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ABSTRACT 


A method  is  presented  for  constructing  finite  difference  schemes  of 
high  accuracy  for  ordinary  differentied.  equations.  Since  it  is  based  on 
power  series  representation  of  the  solution  it  is  applicable  to  partied, 
differential  equations  as  well.  The  method  is  demonstrated  on  the  general 
nonlinear  first  order  ordinary  differential  equation  and  the  truncation 
errors  are  analysed.  A-stability  of  the  method  is  analysed  and  it  is  found 
to  be  stable  in  this  sense  to  arbitrarily  high  order.  Finadly  a nonlinear 
example  is  presented. 
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SIGNIFICANCE  AND  EXPLANATION 


Ihe  author  proposes  a method  for  solving  partial  differential  equations 
numerically.  Hie  present  paper  introduces  the  method  by  applying  it  to 
ordinary  differential  equations.  Ihe  distinctive  feature  of  this  method  is 
that  it  may  be  applied  in  more  and  more  accurate  versions  in  a simple  way, 
and  in  this  respect  it  is  more  easily  used  than  some  classical  methods  such 
as  the  Runge-Kutta  method.  When  applied  to  partial  differentied.  equations 
it  has  no  equived.ent  since  it  is  of  hi^  accuracy  in  both  variables.  Bie 
method  is  based  on  an  old  classical  technique  for  solving  ordinary  differen- 
tial equations  where  power  series  in  x are  substituted  directly  into  the 
equation. 

One  diffictalty  with  numerical  methods  is  that  errors  introduced  in  the 
course  of  generating  solutions  may  grow  so  that  results  may  become  non- 
sensical. When  this  happens  the  method  is  sadd  to  be  unstable.  Hie  pro- 
posed method  is  cuiedysed  with  this  in  mind  and  is  shown  to  be  stable  to  all 
orders  in  a parti culeu:  sense.  Hie  terra  A-stad>le  is  due  to  Dalquist,  and 
indicates  that  a confuted  solution  decreases  when  the  exact  one  does  for  a 
test  equation  that  he  proposed. 


PONER  SERIES  METHODS  I - ORDINARY 
DIFFERENTIAL  EQUATIONS 


Robert  D.  Small 


Introduction 


The  purpose  of  the  series  of  papers  of  which  this  is  the  first  is  to  present  a stable 
finite  difference  scheme  for  peurtlal  differential  equations  that  is  of  high  accuracy  in  all 
variables.  Traditionally,  beginning  with  Euler's  method,  finite  difference  schemes  have 
been  constructed  by  approximating  the  derivatives  that  appe2u:  in  an  equation  by  linear 
combinations  of  function  values  at  equeU.ly  spaced  points  along  a grid.  When  placed  into 
the  equations  the  result  is  a difference  equation  in  the  function  ved.ues  at  grid  points.  A 
major  difficulty  in  taking  a high  order  approximation  for  each  term  in  the  equation  is  that 
the  resulting  difference  equaticai  is  very  often  unstable.  One  way  of  overcoming  this 
problem  is  to  approximate,  to  high  order,  the  solution  rather  than  the  terms  of  the 
equation.  A great  variety  of  methods  known  as  finite  element  methods  fall  in  this  catagory. 
Exanples  are  GcU.erkin  methods  that  fit  solution  surfaces  over  the  cells  into  which  the 
domain  is  partitioned  according  to  some  criterion  that  minimizes  the  error  over  each  cell 
and  collocation  methods  that  minimize  or  annihilate  the  error  at  specific  points  in  the  cells 


These  methods  have  achieved  considerable  success  at  producing  accurate  and  stable  methods 
of  approximating  solutions  to  partied,  differential  equations  but  at  the  ejqpense  of 
conplication  in  their  use  emd  emed.ysis  especially  when  high  accuracy  is  desired  in  ed.1 
variables.  The  author  proposes  eui  alternative  method  that  is  a mild  extension  of  the  power 
series  method  that  is  customarily  used  on  the  ordinary  differential  equations  whose  solu- 
tions are  special  functions.  In  this  respect  it  resembles  Runge-Kutta  schemes.  In  fact 
when  Runge-Kutta  schemes  are  used  for  one  of  the  variables  [1]  in  pareUxilic  equations  it  is 
the  use  of  implicit  schemes  that  promotes  stability  just  as  it  is  in  the  power  series 
method  we  present. 


Sponsored  in  part  by  the  United  States  Army  under  contract  No.  DAAG29-75-C-0024  and  the 
National  Research  Council  of  Canada  under  Gr2mt  No.  A8785. 


Itie  essence  of  the  method  is  to  partition  the  domain  into  rectangular  cells  and  to  fit 
solution  surfaces  over  each  cell  by  means  of  power  series.  Ihe  differential  equation  and 
boundary  conditions  or  continuity  conditions  then  fix  all  the  coefficients  of  the  power 
series,  "nie  point  of  e scansion  of  each  series  is  chosen  in  order  that  the  resulting  scheme 
be  stable  when  the  series  are  truncated  and  this  is  usually  the  center  of  the  cell.  Bie 
accuracy  of  the  scheme  depends  on  the  number  of  terms  of  the  series  that  are  retained  at 
truncation  and  on  the  size  of  the  cells. 

The  purpose  of  this  paper  is  to  present  a stable  form  of  the  power  series  method  for 
ordinary  differential  equations  that  generates  finite  difference  schemes  of  high  order. 

It  has  stability  properties  similar  to  the  implicit  Runge-Kutta  methods  and  accuracy  that 
conpares  favorably  with  those  methods.  Its  advantage,  however,  is  that  it  may  be  straight- 
forwardly applied  to  each  variable  in  partial  differential  equations.  -Bie  present  paper 
lays  the  groundwork  for  the  derivation  of  schemes  and  the  stability  and  error  analysis  in 
partial  differential  equations  that  »d,ll  be  dealt  with  in  a later  paper.  TO  that  end  the 
method  is  applied  to  first  order  nonlinear  ordinary  differential  equations  and  tnsication 

i 

error  is  investigated.  Hjen  the  st2d>llity  of  the  method  in  the  sense  of  Dalquist  is 
studied.  This  detailed  analysis  turns  out  to  be  of  great  iiqportance  when  dealing  with  the 
heat  equation  in  the  paper  to  follow.  Finally  a nonlinear  exemple  is  confuted. 


1.  The  Method 


%/ 


(1) 


Vie  demonstrate  the  method  on  the  first  order  nonlinear  equation 

f(x.y) 


where  f Is  analytic  In  both  x and  y at  some  point  (Xta^).  We  talce  the  solution 


y(x)  of  the  form 

(2) 


® a (x  - x) 
n"0 


nt 


and,  substituting  into  (1)  we  obtain  the  following  hierarchy  of  relations  by  equating  like 
powers  of  x - x. 

a^  - f 

a^  - > a,f‘°'"> 

a3-a2f‘°'l>  tf«'0)  + Ja^f + a2f<°'2) 


(3) 


a^  - a3f<0'l>  . Sa^f'"'^’  + ^a^a^f^^'^)  ^ ,(3,0)  ^ 3^^,(2,1)  ^ 3^2,(1,2)  ^ ^3,(0,3) 


we  have  used  the  notation 


, (n  ,m) 


a"^f(x,Y) 

3x"3y” 


X 

y 


X 

a„ 


This  set  of  equations  evaluates  the  a^  sequentially  end  hence  the  solution  (2)  now  depends 
on  the  unknowns  a^  and  x.  The  analyticlty  of  f is  sufficient  to  define  eOll  the 
from  (3)  2uid  hence  the  function  y (x)  is  analytic  at  x. 

He  may  assuiee  that  the  solution  is  known  at  some  point  x^  and  we  wish  to  find  it  at 
the  point  Xjj  + h.  If  we  choose  x ■ x^  then  a^^  is  known  Inmediately , starting 
seqtmnoe  (3)  then  (2)  may  be  evaluated  at  x^  + h.  This  is  known  as  an  e]q)llclt  method 
since  the  solution  is  calculated  by  successive  substitutions  emd  it  is  precisely  the  power 
series  Mthod  used  to  evaluate  ordinary  differential  equations  expanded  about  analytic 
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points.  As  we  shall  see,  when  (2)  is  truncated  at  some  order  of  accuracy,  this  procedure 
can  be  unstable  in  the  sense  of  Dalquist  [2].  Our  prime  motivation,  however,  is  to  produce 
a method  that  is  stable  when  (1)  is  a peurtial  differential  equation  due  to  the  depenc^nce 
of  f on  further  independent  variables  and  operators  in  these  variables.  In  this  case  the 
explicit  method  is  almost  invariably  unstable  in  the  higher  approximations. 

Tto  eilleviate  the  instedsility  we  place  the  point  of  e]q>euision,  x,  not  at  the  point 
where  the  solution  is  known,  but  elsewhere  in  the  interval  ^ 

and  y (Xq  h)  both  depend  in  a complicated  was  on  a^  which  is  now  a parameter.  Ihe 
solution  yCXq  h)  becomes  known  inplicitly  in  terms  of  y<Xg).  If  we  set 

X = Xq  + Xh 


we  have  the  relations 


= I 


a (-Xh)‘ 


(4) 


0 n! 

n=0 


» a ((1  - X)h)‘ 
y(Xo  + h)  = I 


n=0 

where  the  a^  are  given  by  (3).  Itiese  equations  determine  yCXq  + h)  when  a^  is 
eliminated.  The  io^licit  solution  (4)  may  now  be  approximated  by  truncating  the  series  at 
the  finite  value  n « N,  and  the  solution  is  recovered  by  truncating  (2)  at  the  same 
place , finally  giving  the  set 

N a (-Xh)" 

\ 

n=0 

N a ((1  - X)h)" 
n=0 

N a (x  - x)" 


(5) 


-(x)  = I 


n=0 


n 1 


where  the  a^  satisfy  (3)  for  n = 1,  2,  ...,  N.  Having  produced  a known  value  of  tne 
function  at  x^  + h we  can  repeat  the  procedure  and  determine  the  solution  at  larger  and 
larger  values  of  x by  shifting  all  x values  in  (3)  and  (5)  by  h at  each  step. 


It  happens  that  choosing  x to  be  the  mid-point  of  the  interval  gives  a stable  scheme  for 
the  linear  case  and  ws  analyse  this  in  section  3.  For  the  present,  however,  we  investigate 
the  restrictions  on  f(x,y)  that  lead  to  a workable  schesm,  and  the  acctiracy  that  results 
from  a choice  of  N. 


I 


2.  Truncation  Error  emd  Conaistencv 

It  has  been  noted  that  f(x,y)  must  be  analytic  in  x and  y at  all  points 

found  in  the  course  of  the  confutation  in  order  to  define  the  a . Vfe  consider  the  a 

n n 

as  functions  of  a^  from  (3)  and  regard  a^  as  dependent  on  h and  y(Xjj)  when  (4)  is 
enforced.  Then  to  have  the  a^  bounded  eis  h -►  0 we  require  that  f(x,y)  be  euialytic 
in  X and  y for  all  points  (x^^  + h,  a^)  in  some  range  0 < ii  ^ h.  He  can  now  oosfare 
the  a^  to  the  a^.  By  means  of  (4)  and  (5)  we  have 


(o) 


N a (-Xh)" 

I — 

n nl 
n«0  n-0 


» a^(-Xh)^ 

I — 

« nl 


Since  the  a are  established  by  (3)  we  can  write  a - a (a.)  so  that  (6)  gives  a in 

n n 0 ^0 

terms  of  a^.  He  can  noUoe  an  approximation  to  the  solution  for  small  h by  setting 


*0  + Y (h) 


vdiich  gives 


*n  “ *n‘*0>  Y 


da^ 


1 1 


Placing  these  into  (6)  the  leading  order  solution  for  y is 

,N+1 


(7) 


a^^l(-Xh)‘ 


' (N  + 1) 1 

He  then  obtain  the  result  that  the  a^  differ  from  the  a^  by  terms  of  and  the 

higher  values  of  N introduce  factorials  that  help  out  as  well.  There  is  the  difficulty 
that  (6)  is  a nonlinear  equation  and  the  noticed  solution  for  fixed  h may  not  be  unique. 
It  is  unique  as  h -*•  0,  however,  and  this  may  restrict  the  sise  of  h if  one  cannot  find 
a method  for  identifying  spurious  roots  of  (5) . 

He  now  show  that  the  difference  scheme  (5)  is  consistent  with  the  differential 
equaticm  (1)  as  h -*■  0.  Displaying  the  result  of  substituting  (5)  into  each  side  of  (1)  %<e 
have 


St  , 

N-1 

y 

Vl<* 

dx 

L 

n=0 

nl 

(8) 

f(x,y)  = 

f + 

(fa,o) 

. . • 

+ (Nth  ( 

(a,  f (x,y)  =■  f + (X  - X)  + ... 

~ N 

+ (Nth  ei^zesslon  in  equations  (3))  (x  - x) 

+ higher  order  terns 

nte  arguments  of  f emd  its  derivatives  are  (Xfa^) . Die  "hi^er  order  terms”  are  the 
first  to  feel  the  truncation  of  (5),  the  preoeeding  terms  containing  identical  e3q>res8ions 

to  the  ri^t  sides  of  equations  (3).  Comparing  the  two  e]q>te8sions  in  (8)  the  first 

- , -,N 

SjjCx  - x) 

discrepancy  is  that  the  term  ;;; is  present  in  f (x,y)  but  absent  in  2uid 

NI  dx 


this  gives  an  error  of  0 (h  ) . For  N ^ 1 the  error  vanishes  as  h -►  0 and  we  have  a 
difference  scheme  that  is  consistent  with  the  differential  equation. 

the  error  in  the  solution  is  found  by  coiip2u:ing  the  last  equation  of  (S)  with  (2). 

The  difference  of  the  e]q>ressions  is  of  0(h**^^). 

We  have  foixid  that  f(x,y)  must  be  aneilytic  in  the  range  of  x between  the  point 
where  the  solution  is  known,  x^,  and  the  point  of  expansion,  x,  and  for  all  solution 
values  y that  could  result  there.  Solving  backwards  from  emy  solution  found  for  x to 
the  ri^t  of  the  point  of  e]q>ansian  one  uses  the  same  argument  to  indicate  that  f (x,y) 
should  be  anedytic  for  that  interval  as  well.  When  f(x,y)  has  a point  of  non-^tn^llyticity 
in  X,  that  point  may  be  used  as  zm  end  point  of  an  interval  and  then  the  severity  of  the 
singularity  determines  whether  the  series  (2)  converges.  For  example  if  f(x,y)  changes 
formulas  as  it  does  at  the  Interface  of  regions  with  different  properties,  the  solution  is 
continued  without  difficulty  by  maJclng  the  point  of  interface  separate  two  intervals. 


3.  A-Stability  of  the  Method 


decays  for  all  initial  conditions  when  Re  q < 0.  Haoming  f3]  correctly  criticizes 
indescriair.ant  use  of  this  definition  of  stability  since  errors  cam  still  swan^  the  solution 
they  decay  nore  slowly  than  the  solution  or  solutions  can  doninate  errors  if  they  grow 
faster  than  the  errors.  He  prefer  Dalquist's  definition  over  others,  however,  for  its 
leter  significance  to  partial  differential  equations. 

Since  (9)  is  autonomous  we  can  take  the  solution  of  the  form 

-ax" 

Zn 

-TT- 


and  substitute  into  the  equation  to  obtain 


a . , » qa 
n+1  ^ n 


n • 0,  1,  . . ., 


The  solution  to  this  difference  equation  is  a - q"a  and  hence  we  obtain 

n 


• n n 

y - y a_*_ 

n“0 


idiere  a is  a constant  yet  to  be  determined.  If  the  soluUon  is  known  at  sone  point 
and  we  wish  to  find  it  at  x^j  + h,  then  we  choose  a new  origin  within  this  interval  at  a 
distance  Xh  from  x^j.  Thus  by  a tremslation  of  the  original  co-ordinates  we  may  take 

*0  ” *0  these  values  of  x into  (10)  and  eliminating 

a we  obtain 


He  now  investigate  the  stability  of  the  proposed  method  on  a model  problem  proposed  by 
Dalquist  [21  for  testing  nvmierical  methods.  A method  is  said  to  be  A-stable  if  the 
numerical  solution  to  the  equation 


y(3^  + h)  - Ry(Xg) 


The  first  term  of  (14)  simplifies  further  giving 


(15) 


S / 2N  N 

I + I ? ®nj 

n=0  " „-N+l  j-n-N  n 

for  stability.  We  consider  (15)  as  a condition  on  1 for  given  o,  8 emd  N.  First  we  note 
two  facts  about  the  f^(X),  Jooth  easily  proved  by  induction. 


<0  for  X < - 
2 


2n 


0 for  X«-)n»l,  2,  3, 


>0  for  X > i 
2- 


^2n+l  ^ ° " • 0'  1'  2.  - . 

Using  both  these  facts  we  see  that  the  first  term  of  (15)  is  a sum  of  positive  terms  if 
X > j.  The  second  term  must  be  evaluated  by  brute  force.  We  display  the  first  few  orders  of 
(15)  for  illustration. 


N = 1 2a  + (a^  + 8^)f2  > 0 

N = 2 2a  + 2a^f  - (a^  + oB^)f,  + -^(a^  + 2a262  + B^)f  >0 

A 3 4 4 

N = 3 2a  + - ^a^fj  +^(70“  + Ba^B^  - B^)f^  - i(a®  + 2a^B^  + aB^)f5 

+ + 3a^B^  + 3a^B^  + B®)fg  > 0 

We  recall  that  a > 0 but  g may  have  either  sign.  Then  the  inequalities  for  N * 1,  2> 

are  satisfied  if  X i j.  For  N « 3 the  only  term  that  is  not  positive  for  X i ^ is 

12  scheme  is  to  be  independent  of  a and  B then  we  must  choose  X ” j 

to  suppress  the  offensive  term.  If  the  scheme  may  be  adapted  to  the  given  values  of  a and 
B then  further  stable  values  of  X can  be  found  since  sufficiently  large  values  nuOce  the 
last  term  of  the  inequality  dominate.  These  are  unreasonable  choices > however,  since  X “ 
is  a more  accurate  scheme. 

For  higher  values  of  N we  encounter  more  and  more  terms  whose  signs  wor)t  against  the 

inequality.  As  in  the  case  N = 3,  we  ta)ce  X = j to  suppress  the  terms  with  even 
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subscript  in  f.  It  turns  out  that  for  N ■»  4 this  is  sufficient  for  that  inequality  to 
hold  unconditionally,  but  for  N ^ 5 there  are  always  points  (a, 6)  that  violate  (15).  To 
locate  these  points  we  set  |r|  = 1 and  ^ ^ in  (12)  and  search  for  roots  numerically. 

The  solutions  are  curves  in  the  complex  plane  of  the  variable  qh  that  separate  the  regions 
|r|  >1  and  |r|  <1.  It  was  found  that  writing  qh  in  pol^u:  co-ordinates  gave  the 
simplest  form  for  programming. 

The  results  are  plotted  in  Figure  1.  The  regions  of  instability  are  the  interiors  of 
the  closed  curves  located  in  the  left  half  plane  corresponding  to  the  N of  interest.  The 
several  chains  of  curves  have  been  shifted  vertically  since  they  actually  overlap  at  their 
bases.  To  obtain  an  A-stable  scheme  for  equation  (9)  for  a given  accuracy  N,  we  must 
choose  h in  such  a way  that  qh  does  not  fall  in  one  of  the  undesiredile  regions  when 
Re  q < 0.  This  is  clearly  an  easy  matter  since  the  regions  are  relatively  small  and  sparse. 
Of  course  A-stability,  as  we  have  noted  before,  is  just  a criterion  that  indicates  whetner 


the  computed  solution  decays  when  the  exact  one  does  and  Figure  1 similarly  indicates  that 
if  qh  avoids  the  interiors  of  curves  in  the  right  half  plane,  computed  solutions  will 
grow  only  when  exact  ones  do.  Figure  1 and  equation  (16) , to  be  derived,  will  have  consider- 


able importance  when  dealing  with  partial  differential  equations. 

TO  study  these  cuarves  further  we  rewrite  (12)  in  yet  another  way.  Setting  ^ ^ in 

(12)  we  collect  powers  of  ^ to  obtain 


K 1 y ^ • 

R (2n+l)!  ' 

n— u 


When  |R|  =•  1 the  factor  (R  + 1)/(R  - 1)  is  purely  imaginary  and  the  circle  in  the  R 
plane  maps  into  the  entire  imaginary  axis.  Thus  we  set 


and  ta)ce  in  €uldition. 


to  obtain 


z - f . Ni  = [|]  , = [^J 
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tlwM  curves  srs  the  loci  of  roots  of  the  equartlon 


I 

n«0 


2n 

s 

(2n)l 


2n-i-l 
z 

(2n'«-l)l 


for  reel  psraaeter  c.  Ihe  large  nisdiers  indicate  the  value  of  N for  which  the  curve 
applies,  the  imaginary  axis  is  a locus  for  all  N.  For  cosparison  with  the  text,  z * 

-12- 


1 2n  2 2n+l 

tIktt  “ ifsTITT  • 

n*o  n*o 

The  curves  of  Figure  1 are  the  roots  of  the  Nth  degree  equation  (16)  paraMeterised  by 
reed  c.  Oie  iioaginary  axis  is  also  a locus  of  roots  of  (16)  for  all  N since  any  iaaginary 
z corresponds  to  a reed  c.  Inspection  of  Figure  1 indicates  that  whenever  a curve  has  a 
vertical  tangent  it  intersects  a curve  of  neighboring  order.  This  fact  will  be  proved  fro* 
(16)  and  it  then  serves  to  indicate  that  the  chedns  of  curves  are' infinite  in  extent. 

Since  c paranete rises  the  curve,  a tangent  vector  to  the  curve  is  given  by 

*2  2n+l 

i y -* 

^ (2n+l)  I 


"l-l  2ntl  «2  2n 

y -5 ic  y -£ 

(2n+l)l  (2n)l 


2 2n 

y -5 

L (2n)l 


Ns  write  z • re 


and  introduce  the  notation 


Y r^"co32ne 
i (2n)l 


j Y r^"ain2n6 
l“nio  »n). 


t r^"*^oos(2n+l)e 

I (2n+l)  I 

n“0 


j 2n+l 
n*0 


2 (2n+l)l  ' 2 C 

Then  the  derivative , written  as  real  and  imaginary  peurts  beoosms 


ain(2n-*-l)e 
(2n-fl)  I 


where  [ ] ' 
Written  in  new  notation,  (16)  is 


-^2  1^2  ^°°irSV°2 

, . - ■ 


N N 

Cl  .-cSj 
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n 


•nd  with  c eliminated « the  curve  of  order  N is 


H N N N 

r*  ^ 


Cj  Cj  + Sj  - 0 . 


Similarly  the  curve  of  order  N - 1 is 


and  the  curve  of  order  N + 1 Is 


N N N N 

•=1%  ^ VS  “ ° 


N-  N,  »r  M, 

- s°s'  ■ ° 


where  N3  - and  • 


We  can  specify  points  on  the  curve  of  order  N that  have  vertical  tan9ents  by  having  (19) 
hold  and  by  requiring  the  real  part  of  (17)  to  be  rero.  Eliminating  c from  the  latter  by 


eans  of  (18) , we  )iave 


. c,\»- . 


It  remains  to  show  that  a point  specified  by  (19)  and  (22)  lies  either  on  curve  (20)  or  (21) 
We  consider  separately  the  cases  of  N even  and  odd.  For  N - 25,  we  have  - 5, 

Wj^"N,  H2“N-1,  Hj^H  - 1.  Then  (19)  end  (22)  can  be  written 

Ur<f7w/l/’ 

Either  *0  or  the  trivial  solution  c”  - 0,  - 0 holds.  The  former 

is  (20)  and  the  latter  satisfies  (21).  Then  for  N - 2N  -f  1,  we  have  - 5 + 1,  ■ 5, 


M,  ■ N,  M « N - 1,  emd  (19)  and  (22)  become 


< <\K\  « 


-r‘  -rv  w 


N 1^1  M N N 

Either  “0  or  C^  ■ 0,  S^  ■ 0.  Ihe  fomer  is  (20)  and  the  latter  satlsflea 

(21). 

The  preoeeding  gives  some  evidence  that  the  chains  of  curves  of  Figure  1 extend 
indefinitely.  It  Is  iaposslble,  for  exasfjle,  to  have  an  Isolated  siaooth  closed  curve  that 
could  be  Blssed  in  the  nmerical  search  for  such  curves,  since  any  point  where  the  tangent 
is  vertical  must  necjessarily  lie  <xt  another  curve  m well. 


4.  A tionllnear  Exanple 


The  equation 


uMd  by  Rosser  [4],  is  chosen  to  illustrate  the  accuracy  of  the  method.  In  anticipation  of 
treating  partial  differential  equations  by  this  nmthod  we  used  the  inc>licit  scheme  with 
1 B ^ since  this  is  the  most  stable  form  of  the  difference  scheme.  For  strict  economy  of 
effort  in  solving  (23)  the  ejqtlicit  scheme  with  X » 0 would  be  better  since  it  is  stable 
for  small  enough  h in  the  sense  that  errors  do  not  dominate  the  solution  and  it  is  solved 
by  substitution.  The  coefficients  (3)  were  derived  for  equation  (23)  for  N • 1 to  5 and 
system  (5)  set  vp.  The  first  equation  of  (5) , 

(24)  I 2 

n-0 

is  a nonlinear  equation  for  a^  and  must  be  solved  by  iteration.  Newton's  amt)x>d  was  used 
since  the  derivatives  of  a^  with  respect  to  a^  are  easily  fowid  from  formulas  (3). 

Since  for  large  values  of  h the  solution  can  <Aiange  a great  deal  in  tJie  distance  ^ 
between  the  point  where  the  solution  is  Jcnown  md  the  point  of  mjqpansion,  the  first 
guess  for  Newton's  method  really  should  not  be  t)ie  iMt  )cnown  solution  value  but  the  much 


? 


^computed  ^exact 
^exact 

was  tabulated  <uid  usually  occurred  at  x ~ 0.  There  are  present  in  the  tedile  some  random 
effects  and  some  effects  due  to  symmetry  about  x ~ 0.  Generally,  however,  one  can  detect 
the  improvement  in  accuracy  as  h decreases  and  as  N increases. 


h 

N 

- Order 

of  Method 

« 

Step-size 

1 

2 

3 

4 

5 

2 

OP 

-6.5 

X 

1 

o 

OP 

2.1 

-2.9 

X 10“^ 

1 

OP 

-4.4 

X 

lO"^ 

4.3  X 10"jj 

-2.6  X lO"^ 

1.2 

X 10“^ 

.5 

00 

-2.0 

X 

lo”^ 

1.3  X 10-^5, 

-1.5  X lO"^ 

-5.8 

* “usi 

.25 

2.7  X lo"^ 

-6.3 

X 

O 

1 

-7.9  X 10“® 

-7.5  X lO"® 

1.8 

X lO"^ 

.1 

3.3  X lo”^ 

-1.1 

X 

1 

o 

-6.5  X lo”® 

-1.8  X lO”* 

-3.6 

* ^°(.3) 

Table  of  Maximum  Relative  Errors 


The  nuucimum  relative  errors  over  all  grid  points  are  tabulated  for  the  solution  of 
2 1 

-2xy  , y(-8)  » — in  the  range  (-8,«>)  for  various  step  sizes  and  orders  of  the  method, 
od 

error  of  <»  indicates  that  at  some  step  the  difference  scheme  failed  to  have  a real  solu- 
tion. If  the  maximum  error  did  not  occur  at  x = 0,  that  grid  point  appears  in  brackets. 


Concluding  Reaarks 

In  suawtry,  the  povrer  series  method  is  offered  as  a way  of  discretizing  partial  dif- 
ferential equations  with  high  order  of  accuracy  in  all  independent  variables.  This  intro 
ductory  paper  demonstrates  the  discretization  procedure  on  first  order  nonlinear  ordinary 
differential  equations  and  analyses  A-stability  of  the  method.  An  implicit  form  of  the 
method  turns  out  to  be  A-stable  to  all  orders  subject  to  a mild  restriction  on  the  step- 
size.  The  key  equation  <16)  that  determines  that  restriction  arises  frequently  when  deal 
ing  with  simple  partial  differential  equations.  An  example  is  the  heat  equation  that  is 
dealt  with  in  the  following  paper. 
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